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Abstract 

In a recent letter [Phys. Rev. Lett. 100, 164101 (2008)] and within the context of quantized 
chaotic billiards, random plane wave and semiclassical theoretical approaches were applied to an 
example of a relatively new class of statistical measures, i.e. measures involving both complete 
spatial integration and energy summation as essential ingredients. A quintessential example comes 
from the desire to understand the short-range approximation to the first order ground state contri- 
bution of the residual Coulomb interaction. Billiards, fully chaotic or otherwise, provide an ideal 
class of systems on which to focus as they have proven to be successful in modeling the single 
particle properties of a Landau-Fermi liquid in typical mesoscopic systems, i.e. closed or nearly 
closed quantum dots. It happens that both theoretical approaches give fully consistent results 
for measure averages, but that somewhat surprisingly for fully chaotic systems the semiclassical 
theory gives a much improved approximation for the fluctuations. Comparison of the theories 
highlights a couple of key shortcomings inherent in the random plane wave approach. This paper 
contains a complete account of the theoretical approaches, elucidates the two shortcomings of the 
oft-relied-upon random plane wave approach, and treats non-fully chaotic systems as well. 

PACS numbers: 03.65. Sq, 05.45.Mt, 71.10.Ay, 73.21.La, 03.75.Ss 



1 



I. INTRODUCTION 



Finite and low-dimensional quantum systems often possess statistical properties whose 
deviations from universality contain some basic dynamical information about the system 

. A recurring challenge is to understand precisely what information is buried in those 
statistical deviations and how to extract it. Before that can be addressed, the universal 
behaviors themselves must be understood. At the heart of many such universalities generally 
lurks a connection to the Bohigas-Giannoni-Schmit (BGS) conjecture , which asserts 
that systems with underlying chaotic dynamics have the fluctuation properties found in 
random matrix theory [3] • If interest lies in quantities for which the position representation of 
the eigenf unctions is critical, a random plane wave model is often introduced that augments 
the BGS conjecture js], ol. The primary goal of extracting system specific information can 
then proceed, but generally requires a more powerful theory. 

The preponderance of statistical measures heretofore introduced for analysis are local 
in energy, configuration space, or both. Examples are given by the Dyson-Mehta cluster 
functions [?!, and the amplitude distribution and short-range two-point correlation function 
c(|r — r'l) = {ilj{r)ilj{r')) of a given eigenfunction, '?/'(r). On the other hand, there has 
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been a rather recent introduction of new, non-local statistical measures 1^, 
have been motivated by the need to understand the interplay between interferences and 
interactions in mesoscopic systems. For example, one place where this interplay is known 
to have an important role is the addition spectra of quantum dots; other examples are 
coming from cold fermionic gasses. Focussing on just one motivation, the addition spectrum 
is experimentally accessible throu gh the position of the conductance peaks in Coulomb 



Blockade transport measurements [l3|, |l4 
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161 . llTj . Fluctuations of peak spacings are 



associated with interference effects, and experimentally are clearly incompatible with a non- 
interacting description of the conduction electrons in the dots. Indeed, a non-interacting 
description predicts a strong bimodality of the peak spacing distribution - associated with 
an odd-even character of the number of (spin-1/2) electrons - and this is not observed. 
Assuming further that the dot possesses a chaotic dynamics, the distribution for odd-A^ 
spacings should show a characteristic Wigner surmise shape, whereas a density similar to a 
Gaussian with extended tails is observed. 

At first, it was argued that a no n- Fermi liquid description of the interacting electrons 
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might be necessary 



now emerged 
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;o in terp ret the experimental data. However, the picture which has 
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221 ] is that although an understanding of both peak spacings 
and ground state spin distributions is still incomplete, it is reasonable to expect that most 
phenomena will eventually be explained within a Fermi liquid framework. More specifically, 
the electrons can be thought of as quasi-particles confined by a potential f/conf(r), which 
could in practice be computed within a self-consistent Thomas-Fermi-like approximation 



231]. They also interact weakly through a screened Coulomb interaction Vsc{r,r'). For dots 
significantly larger than the screening length, confinement does not modify appreciably the 
screening process, and the bulk expression for Vsc{r — r') is appropriate. Furthermore, for 
the experimentally relevant gas parameter being of order one, the screening length is not 
much different from the Fermi wavelength A^.. Under these circumstances, the problem is 
well described by the short range approximation 

V;,(r-r') = ^5(r-r'), (1) 
with 1/ the mean local density of states (including the spin degeneracy) {v = m/Tifi? for d = 2) 



and Fq the dimensionless Fermi liquid parameter [2J] (for d = 2 and of order one, F^ is 
in the range 0.6-0.8). Boundaries could potentially modify this picture somewhat, but it is 
expected that a slightly modified Fq would be sufficient to capture the effects; since this is 
not the focus of our study, we leave it for future consideration. 

In this approximation, the first order contribution of the residual interactions to the 
ground state energy can be expressed in the simple form 

SE''' = ^ J rfrnT(r)n;(r), (2) 

with no- the unperturbed ground state density of particles with spin a. This expression 
was the starting point for a study demonstrating the increased importance of 6E^^ if the 
dynamics are not fully chaotic jlO]. 

Because n^{r) can be expressed as a sum over the absolute square of single particle 
eigenfunctions which are occupied, the mesoscopic fluctuations of the residual energy term 
6E^^ of Eq. ([2]), or of similar quantities, are related to the fluctuations of the one-particle 
eigenstates of the unperturbed system. However, in contrast to the correlation function 
c(]r — r']), Eq. ([2]) involves both an integration over space and summation over energy, and 
therefore in this way, it is probing new aspects of the fluctuation properties of the eigenstates. 



The goal of this paper is to follow up on our recent work [25] on the average and fluctuating 
parts of the quantity in Eq. ([2]) in several ways. To begin with a complete account is 
given of two theoretical approaches using random plane waves and semiclassical theory. As 
already shown, for fully chaotic quantized billiards the two methods give identical leading 
functional dependence on wave vector and system size for average and fluctuation properties. 
In addition, the significant differences between Dirichlet and Neumann boundary conditions, 
can be understood. On the other hand, in the statistical limit of uniformity, the fluctuation 
properties differ in two ways in the pref actor. The semiclassical treatment in the spirit of the 
Gutzwiller trace formula [2^, helps identify dynamical correlations and a term missing 
from the expressions derived with the random plane wave model. Thus, the Gutzwiller 
periodic orbit approach provides both a deeper understanding of the mechanism underlying 
the fluctuations and a good quantitative agreement with exact numerical calculations for 
the examples considered. Next, non-fully chaotic systems are partially addressed. The 
fluctuations are found to be greatly magnified in two essential ways. One relates to non- 
uniform projected classical densities and the other to enhanced dependence on wave vector 
and system size. 

The organization of the paper is as follows. In Sect. II, the necessary background material 
and notations are introduced, including a more precise definition of the statistical quantities 
to be studied. Also introduced in this section is the random plane wave model, which is then 
used to analyze in Sect. Ill the mean and fluctuating behaviors. The random plane wave 
modeling reproduces quite accurately the mean behaviors, but it fails by most of an order 
of magnitude to predict quantitatively the fluctuations. This motivates the introduction in 
Sect. IV of the semiclassical analysis in terms of classical trajectories. This approach corrects 
the random plane wave method overestimate of the fluctuations. In Sect. V, extended 
semiclassical methods are introduced for systems which are not fully chaotic. Finally Sect. VI 
contains a discussion and summary. 



II. PRELIMINARY CONSIDERATIONS 



To begin, it is worth motivating the introduction of a non-local statistical measure in a 
little more detail. For example, consider two-degree-of-freedom chaotic quantized billiard 
systems with one body eigenstates ipi{r) and energies Ei. In the absence of interactions, the 
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many-body eigenstates are Slater determinants characterized by spin-dependent occupation 
numbers fi^a = or 1. Within the short range approximation, Eq. ([1]), the contribution 
of the residual interactions to the ground state can be written in first order perturbation 
theory as 

^^Ri = ^J2 /.(+)/.,(-)^^. (3) 

where A is the mean single particle level spacing in the neighborhood of the Fermi surface. 
Mij is given by 

M,,=A J rfr|^,(r)H^,(r)|2 , (4) 

where A is the area of the billiard. With this definition, the {Mij} are dimensionless 
quantities with a mean value expected to be roughly equal to unity for i ^ j, and to three 
for i = j. These expectations would apply to uncorrelated Gaussian random amplitudes 
assuming time reversal invariance holds; ahead more precise results are derived. 

Of main interest is how the value of 6E^^ changes when a particle is promoted from one 
orbital to another (as opposed to the variations of the residual energy as particles are added 
into the system). Consider that the ground state of the non-interacting A^-particle system is 
such that the levels below the Fermi energy are doubly occupied except for the last level ip, 
which may be singly or doubly occupied depending on the parity of A^. Promoting a particle 
from the orbital ip to ip + 1 has a one-particle energy cost (-Ej^+i — Ei^). If however this is 
compensated by the corresponding difference in residual energy, the ground state occupation 
numbers /j,(±) will be modified by the interactions, yielding in some circumstances non 
trivial, i.e. different from or 1/2, ground state spins. Imagining {/j (_|_)Mjj/j (_)} in the 
form of a (square or nearly square) symmetric matrix, shifting an occupancy from one level 
to another subtracts the column or row being vacated and adds a column or row to the 
newly occupied orbital (row or column depends on the spins of the removed and added 
particles). Apart from a couple individual Mj^'s near the diagonal, the difference in the 
residual interaction can therefore be expressed in terms of the difference of two sums of the 
form 

i 

S, = J2 M., . (5) 

As the {Mij} are positive definite. Si has a locally defined, increasing- with- index, positive 
mean. However, because one column or row is added and another subtracted, their means 
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largely cancel, and thus the mean Si behavior cannot be involved in altering the ground state 
occupancies of the single particle levels defined by f/conf(i")- [The mean of the few individual 
Mjj's near the diagonal not included in Si might though.] However, if the fiuctuations of Si 
or {Mij} are sufficiently large, they have the potential to alter the nature of the ground state. 
Thus, statistical measures based on the properties of the Si are of fundamental interest, in 
particular their mean values and fiuctuations. 

The Si have the unusual character that they are integrated over space and involve a 
sum over eigenenergies. Their investigation thus requires two basic ingredients. With the 
definition 

Nir;E)= dE'nir; E') = Tl^iir)^ 6 (E - Ei) (6) 
•^0 i=l 

Si can be expressed as 

S, = Ajdr \^,ir)fJ2\'^,ir)f = A J dr |v[/,(r)|' iV(r; , (7) 

j<i 

with the understanding that Ei < Ef < -Ei+i (and assuming for simplicity that there are 
no degeneracies). One of the two required ingredients is the behavior of N{r;E). As it 
results from a summation over the absolute square of eigenstates up to a certain energy, 
it is dominated by a secular behavior; see Fig. 13 in Ref. for example. The secular 
component Nsec{T^] E) emerges from an energy smoothing which, although local, is also 



necessarily broader than the Thouless energy 291]; ahead the notation (■) is introduced 
to denote this averaging. This energy smoothing implies that only dynamics on a time 
scale shorter than the shortest periodic orbit is relevant, and thus this decomposition is 
independent of whether the system dynamics is regular, fully chaotic or has some other 
character. A'^spnfr; E) has been shown to be given by an excellent semiclassical (asymptotic) 



approximation 
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N^ec{r;E} = 



kx 



(8) 



A 

where the coordinate x is defined locally as the perpendicular distance from the boundary, k 
is the magnitude of the wave vector at energy E, and the + sign is for Neumann boundary 
conditions and the — sign for Dirichlet. Here N\y{E) refers to just the leading term of the 
Weyl formula, N\y(E) = ^^E. The validity is governed by A;L ^ 1, where L is a length 
scale, specified ahead in the paper, but which necessarily must be shorter than the width of 
the system. Away from the boundary, the secular behavior approaches an overall constant. 



However, the existence of boundary conditions and a minimum wavelength scale combine 
to create persistent oscillations (Friedel oscillations), which are maximal near the boundary 
and which fade away toward the interior of the billiard. The negative sign of the Bessel 
function for Dirichlet boundary conditions respects the vanishing of eigenfunctions at the 
boundary as it must. Note that just as the critical portion of the density of states can be 
expressed as a series with volume, boundary, curvature, and oscillatory components, the 
same is true of A^(r; E). The above expression does not include the curvature components 
and must therefore be missing at least part of the 0{[kL]^^) corrections. 

As noted Ai'sec(r; -E) is a smooth function of the parameter E, but Si actually involves 
A^sec(r; E^), which changes abruptly at the points where the energy surpasses each eigenvalue 
(i.e. is a function of i). The former quantity does not have a monotonous dependence in the 
number of particles i since Ei contains the Gutzwiller corrections from periodic orbit theory 



26l . 127(1 that determine the precise positions of the levels. We therefore consider instead the 



slightly modified and properly normalized decomposition 

iV(r; K+) = iV,,,(r; E+) + 5N{t; Ef) 

Ji{2kir) 



N,ecir;El) ~ 



1 ± 



(9) 



where the C is the billiard perimeter and not to be confused with length scale L mentioned 
above. This decomposition has the further advantage that the fluctuations 6N{r; Ef) not 
contained in the secular behavior average to zero when integrated over space. In this way, 
density of states oscillations, which are not of interest here, do not get intertwined with the 
fluctuations that are the focus of this study. To the order of corrections incorporated in 



Eq. ([9]) , ki can equally be defined as \/2mEi //i or as the mean value obtained from the Weyl 
series. 

The other main necessary ingredient is the behavior of |\E'.j(r)|^ which leads to the two 
principal approaches contained in this paper. One approach is to rely upon a statistical 
model which uses an ensemble of random plane waves to mimic the properties of chaotic 



Bogomolny 



eigenstates[8|j9|, ImI] and the other is to use a semiclassical theory building on the work of 



We begin with the random plane wave modeling as it is technically simpler. 
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III. RANDOM PLANE WAVE MODELING 



m 



The random plane wave model [8|, |9|, |31| has been introduced in which the eigenstates 



are represented, in the absence of any symmetry by a random superposition of plane waves 
J2i O'l exp(zki -r) with wave- vectors of fixed modulus \ki\ = kp distributed isotropically. Time 
reversal invariance may be introduced as a correlation between time reversed plane waves 
such that the eigenfunctions are real. Similarly, the presence of a planar boundary imposes 
a constraint between the coefficients of plane waves related by a sign change of the normal 
component of the wave- vector k;. Near a boundary, and using a system of coordinates 
r = X + y with x and y the vectors respectively perpendicular and parallel to the boundary 
(of norm x and norm y), eigenfunctions are therefore locally mimicked statistically by a 
superposition, 

i>i{r) = —— aics {ki ■ x) cos (k/ ■ y + cpi) (10) 

™ 1=1 



where cs(-) =^ sin(-) for Dirichlet and cos(-) for Neumann boundary conditions 33|. The 
phase angle (pi, the real amplitude ai, and the orientation of the wave vector k/, are all chosen 
randomly. The amplitudes a; are zero-centered independent Gaussian random variables with 

{aiai') = dwa"^. 

To complete the model, it is necessary to determine the variance which is fixed by the 
normalization of the wavefunctions. Here this constraint is imposed only on average, rather 
than for each individual state. A priori, proceeding in this way might be expected to miss 
weak correlations between the eigenfunctions. The question is whether one should expect 
them to be insignificant. In principle, the answer is yes, but only if local properties of the 
eigenfunctions are being considered and the effective dimensionality is large as it would be 
for kpL ^ 1, where fci? is the Fermi wave vector. This issue is further discussed in the 
semiclassical theory of Sect. HVl ahead. 

Using the random plane wave representation Eq. ffTOj) we have 

l^i(r) 1^ = --2- ^ aiavcs (k; ■ x) cs (k// • x) cos (k/ • y + (fi) cos (k// ■ y + (fv) , (11) 



efT 



Ll'=l 
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whose expectation value is given by 



a 



eff 1=1 



^ [1 ± COS (2k, • x)] 



47Veff 



1 ± 



27r 



2t7 



d9 cos (2/cjrX cos ^) 



4iVeff 



[1 ± Jo(2A;^x)] , 



(12) 



where we follow the convention that the upper sign refers to Neumann and the lower sign 
to Dirichlet boundary conditions respectively. The ensemble transformation — > 
^ d9i has been employed to simplify the calculation. Above, the norms of the wave 
vectors, equal to /cj, are assumed at or near enough the Fermi surface, that they can be 
denoted by kp. Integrating over the area of the billiard to fix the normalization gives 

/ 2 r rco 
dr {\Ur)f) = ^ A±£J drJo{2kpx) 



4Areff 



1 ± 



2kpA 



(13) 



which fixes the variance o"^ to next to leading order in kpL. 



A. Average Properties 

The first step in calculating the average behavior of Si is to isolate its secular and fluc- 
tuating behavior. The model above implies 

A{M^)\') = -( ^-—^[l±J,{2kpx)] . (14) 

and that is consistent with A^sec(i'; E^), i-e. the Friedel oscillation contributions to A'"(r; E^). 
This form apphes more generally than the random plane wave model. Just as A^sec(r; E^) is 
independent of system dynamics, so also is this result for the same reasons. For example, it 
would emerge for integrable systems as well assuming the averaging is over energy intervals 
greater than the Thouless energy. Finally, note that from the correction to unity of the 
leading constant in this expression, one sees that if the boundary conditions arc Dirichlet, 
the local mean behavior of (|'0i(r)|^) well into the interior is slightly elevated above 1/^ to 
compensate for the reduced density near the boundary, and just the opposite for Neumann 
boundary conditions. 
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Let 



e,{r) = ^|^,(r)r 
1 



1 ± 



2kpA 



-A{\Ur)f)-l 
Jo{2kpx) 



(15) 



Under spatial integration, ej(r) as well as both halves of the second expression above each 
separately vanish whereas only the second half of the expression is affected by taking the 
expectation value and in that case it vanishes. Ahead, this decomposition simplifies the 
discussion of the fluctuations. 

Returning to the calculation of iSj, substituting the relations from Eqs. (PfTSl) . integrating 
the constant terms, and dropping second order terms gives 



5,. 



i ± 



A 



and therefore 



{S^) 



i± 



A 



Ji{2kpx) 
dr — ei(r) 

rv pX 



Ji{2kpX) . , 
dr — ^ (ei(r)) 



(16) 



'^A 



dr 
2£ 



Ji(2kpx) 



Jo{2kpx) 



irkpA 



(17) 



where we choose kp = a/ 2mEi / h. Note the first correction is independent of whether 
the boundary conditions are Neumann or Dirichlet (see Appendix |X] for the calculation of 
second order terms not related to curvature and discontinuities in the boundary). We stress 
furthermore that Eq. f|T7|) is applicable independently of the nature of the dynamics, and 
in particular apply equally well to integrable and chaotic systems. A simple semiclassical 



proof of this will be given in section 
The well known chaotic cardioid 



nYi 



34 



35 



36| and stadium billiards [371, l38|, l39[ are highly 



suited to illustrating the precision of this relation. The two symmetry-reduced billiard 
boundaries are illustrated in Fig. [1] along with their shortest periodic orbits to which we 
return ahead in the discussion of the Fourier transform of the fluctuations. Comparison to 
Eq. ( fT7|) is shown with a computation of Si using the first 2000 odd parity eigenstates of 
the cardioid billiard and the same number of even-even eigenstates of the stadium billiard; 
the latter calculation tests the effects of Neumann boundary conditions. Figure ([2]) plots 
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(a) 




(b) 




FIG. 1: Drawing of the desymmetrized cardioid and stadium billiard boundaries. The dashed lines 
illustrate the paths of the shortest periodic orbits for each system. Whereas this orbit is isolated in 
the cardioid billiard, in the stadium billiard it is a member of a continuous one-parameter family 
of identical orbits, indicated by the grey-shaded rectangular region. 

the differences, Si — (Si), versus the state index for both billiards. As often happens with 
semiclassical approximations, even though the result is asymptotic, it is valid right down 
to either ground state. Note that the second order corrections have not been included in 
the secular behavior equations, Eqs. ([9], [Hj), and so it is seen that the cardioid results are 
not centered on zero, but on a constant somewhere nearby. The same is also true for the 
stadium, except that the mean constant was subtracted in order to compare with the solid 
line predictions from Sect. IV Bl ahead. 

B. Fluctuations 

Consider now the fluctuations of Si. The quantity that actually sets the scale for the 
fluctuations in the residual interaction energy is, as discussed in the introduction, approxi- 
mately the variance, Var[5j — iSj]. It is understood that do not differ by more than some 
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FIG. 2: Comparison of Si and the approximation given in the last hne of Eq. ()17p . The difference 
Si — {Si) is plotted versus the eigenstate number i. In (a), the results for the odd parity eigenstates of 
the cardioid billiard using Dirichlet boundary conditions are shown. In (b), the same results for the 
even-even eigenstates of the stadium billiard with Dirichlet boundary conditions are shown, except 
that the mean overall constant term has been numerically subtracted. Note the significant increase 
in the scale of the fluctuations about the mean for the stadium billiard. Ahead in Sect. IV Bl a 
rudimentary theory is given for bouncing ball modes that leads to Eqs. ()75|76p shown as the dashed 
and solid lines respectively in (b). 



small integer. More specifically, the interest is in computing this quantity to the leading 
order in the semiclassical parameter (kpL). 

Using Eq. ([7]) along with the decomposition into a secular and fluctuating part of A^(r, E^) 
given by Eq. IQ, results in Si being written as the sum of two independent terms. The second 
one, been considered in [l9|, |2l| within the random plane wave 

approximation. It has a variance scaling as \og{kpL)/{kpL), which is of lower order than 
the leading behavior of Var[5j] calculated ahead. Although, the second term is potentially 
of physical interest (see Sect. IVl|) . the focus here is on developing the theory that gets the 
leading term analytically. Applying Eq. (fT5|) and dropping all the lower corrections gives for 
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the covariance between Si and Sj 

Covar(5.5,) = -||dr,|dr, ^'f^'^^'^ ^'f^'^f [{e.{r,)e,ir,)) - (e.(r,)) (e,(r,))] . 

(18) 

As before, x and y are the coordinates perpendicular and parallel to the boundary. Not 
surprisingly, since eigenstate-to-eigenstate correlations are not included in random plane 
wave modeling, 

(e.(ri)e,(r2)) - (e.(ri)) {e,{r,)) = 6,, [(e.(ri)e.(r2)) - (e.(rO) (^(r^))] (19) 

and thus, Var[5j — Sj] = 2Var[5j]. 

Performing a little more algebra gives 

(e,(r0e.(r2)) - (e,(rO) {u{r,)) = A' [([^.(rOn^.lr^)^ - (|^.(rOO ([^.(r^)^] (20) 

where Eq. f fTOj) is applied to evaluate the right hand side of this equation. Each resulting 
term has a product of four Gaussian random coefficients. The fluctuations are thus given 
by pair-wise correlating coefficients such that 

{aiai> ) = {aiai'){amam') + {aiO'm){ai'am') + {aiam'){ai'am) 

= Cr^ (Su'Smm' + SlmSl'm' + Slm'Sl'm) (21) 

where [1,1') are linked to the first coordinate, ri, and (m, m') are linked to the second 
coordinate, T2. The first term, which correlates the wavefunctions taken at the same position 
just reproduces the mean (|?/'j(ri)p) (|'?/'j(r2)P) and cancels from Eq. fl20|) . The two remaining 
terms give the same contribution, which can be understood as a consequence of time reversal 
invariance. Only one of those terms would be non-zero for a time reversal non-invariant 
system, and the result for the variance of Si would just be divided by two in that case. 
Therefore, together with averaging over {(fi,(f2), the variance is 

^" l,m=l 

xcs(k„ ■ xi)cs(k.m ■ X2) cos [k/ ■ (yi - 72)] cos [k^ ■ (yi - 72)] . (22) 

Refiection of either of the vectors (ki,'km) leaves the integrand unchanged. Thus, 
cos [ki ■ (yi - ys)] cos [k^ ■ (yi - ya)] is equivalent to exp [i(ki - k„) ■ (yi - ys)] and can be 
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replaced in the integrand. The summations can be replaced by angular integration again as 
was done in Eq. f|T2|) . 

For the purpose of understanding the asymptotic limit, one is tempted to extend the limits 
of integration in these integrals. However, that generates divergences associated with large 

= Hi ~ y2 and small sin{6i — 6^) (where m is the angle between the vector and the 
direction ct) . This indicates that over large distances the random plane wave model as given 
by Eq. ffTOl) cannot be applied. One way to think of this is to imagine a true eigenstate of 
some chaotic billiard. Locally, one could project onto the form of Eq. (fTOj) and approximately 
solve for a set of coefficients {«/} and plane wave orientations. However, the solution set {a/} 
would be dependent on the location along the boundary where the projection was performed 
due to the rotating orientation of the local coordinate system. Even if the Gaussian random 
modeling were perfectly fine from state-to-state, as r2 got further from ri, the two cross- 
terms in Eq. (!2T|) that generate the variance would progressively decay on a length scale 
given by the typical dimension L of the system. This is related to the behavior observed 
for the spatial autocorrelation function; see Fig. 2 of [40]. (Note that the first term, which 
reproduces the square of the mean would on the other hand not decay.) Ahead, it is seen 
that the results depend only logarithmically on this parameter for the Neumann boundary 
conditions and not at all for Dirichlet, so that it is not necessary to describe very precisely 
this decay as long as the proper length scale is introduced. 

A Gaussian form exp(— 5?/^/2L^) is convenient and gives 

Var = ^ f dOi f dOm f dx '^^ ^ ^^^'^ cs{kpX cos 6i)cs{kpX cos 6„ 
n^A^ y_^/2 J-n/2 Uo kpX 



X/: y d{by) exp (^-^ + «(ki - k„) ■ (yi - y2) J . (23) 

Including the Gaussian cutoff, and noting that the dominant contributions come from regions 
in which sin 5^ is small (with 59 = 9i — 6m), it is possible to approximate sin^^ — sin^^ ~ 
cos{6)66 with 6 = {6i + 6m) /"i. The integrand I{d, 69) becomes 



1(6, 66) = V2ttLC 



2 



l±|sin(^)| ' {kpL cos 969)' 



exp 



(24) 



2/u^ 

where the sign — and + correspond to Dirichlet and Neumann boundary conditions respec- 
tively. Performing the integration over the variables 69 and 9 yields 

Var[5.] = ^{m)e (25) 
14 



where we have introduced the function 




(26) 



[ d{sm9)\\9) = (21n2-l) 



Dirichlet 



(27) 



= (2 In 2 - 1) +4 In 



TikpA 



Neumann . 



(28) 



2£ 



Note caution must be exercised in evaluating the Neumann case. There the angle 6 cannot 
be allowed to decrease to less than the inverse of kpL where the cutoff expression becomes 
invalid. In the absence of other considerations, a very reasonable choice for L is just half 
the average length between two reflections; see Appendix [Bl For a 2D concave billiard, this 
gives exactly L = 'kA/{2C), and this value has been substituted into the Neumann form. 

Finally, consider the difference between Dirichlet and Neumann boundary conditions. 
Equations fl27f28p both roughly imply a kpC behavior. The prefactor ln(2/ Y^)/27r'^ ^ 0.0031 
is rather small in the first case, whereas for Neumann boundary conditions there is a logarith- 
mic enhancement which can be understood as a (much larger) effective prefactor (a factor 40 
larger for i = 1000). From the point of view of the calculation, the difference between these 
two cases can be related to the sign change between 1 — | sin(^)| and 1 + | sin(^)| in Eq. (12^ . 
in such a way that whispering gallery modes (for which the corresponding classical orbits 
have 6 ~ 7r/2) are suppressed for Dirichlet boundary conditions whereas they dominate 
(because they are less affected by the exp[— {kpL cos 66 2] factor) in the Neumann case. 
This makes sense since the main source of Si fluctuations originates from the wavefunction 
fluctuations' probability density |^/'j(r)p in the mean field oc [(1 ± Ji{2kpx)) / {kpx)] gener- 
ated by the Friedel oscillations of all the other particles below the Fermi energy. Dirichlet 
boundary conditions however impose that |'?/'j(r)p ^ as r approaches the boundary and 
therefore inhibits this contribution. 

Figure [3] illustrates the comparison between the analytical results of the random plane 
wave model for the variance, Eq. (1251) . for the cardioid and quarter stadium billiards. For 
the stadium, kpC is replaced by kpC^, i.e. the length of the straight edges where Neumann 
boundary conditions are imposed (even-even symmetry class). The theory for the Dirichlet 
case, cardioid billiard, appears to be roughly a factor six too great. In order to understand 
the discrepancy, the more powerful approach of semiclassical theory is developed in the next 
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FIG. 3: Variance of 5j for the cardioid (odd parity only) and quarter stadium (even-even symmetry 
only) billiards using Dirichlet boundary conditions. In (a), the discrete points are the cardioid 
bilhard results, the dotted line is the result of the random plane wave model, i.e. Eq. (j25p . and 
the long dashed line is the semiclassical theory, Eq. (j52l53p . given in Sect. HVT] ahead. In (b), the 
discrete points are the stadium billiard results, the dotted line is the result of the random plane 
wave model, Eq. ()25p . and the long dashed line is the prediction Eq. (j77p from the semi-quantitative 
semiclassical theory developed in Sect. IV Bl ahead for the bouncing ball modes. 



section. For the stadium, both Neumann boundary conditions and bouncing ball modes 
must be considered. This involves additional complications treated in Sect. IV Bl ahead. 



IV. SEMICLASSICAL APPROACH 



It is important to develop a semiclassical approach. It gives a more powerful theory and 
sheds some light on the difficulties that the random plane wave model is having in providing 
a quantitative description of the Si fluctuations. The most immediate conceptual difficulty 
in getting started is that a treatment of the Si implies, through Eqs. (I15fl6l) . addressing 



16 



the fluctuations of individual wavefunctions, wliereas semiclassical approximations valid for 
chaotic systems - such as the ones based on the semiclassical Green functions - converge 
only for quantities smoothed on an energy range containing a significant number of levels. 
Here, however, this difficulty can be overcome. 



For this purpose, let us, following Bogomolny j32|, introduce a local energy averaging 
that is generally much narrower than the Thouless energy 

= ^ E 5. (29) 

E-^<E,<E+^ 

The notation AA^ =^ N{E + ^) — N{E — ^) represents the number of levels in the energy 
interval [E — E + ^]. With this notation, the variance is 

Var[ San] = ({San- (S))') = -r^Var[5.] + ^^-llcovar,^, [5.5,] (30) 



AN ' ' AN 

((iS) =^ (Si) = ( San))- Computing the locally averaged quantity San, for which conver- 
gent semiclassical approximations can be used, it is possible to extract the variance and the 
covariance of the Si from the scaling in AA^ of Var[ iSaat]- In addition, as expected from 
the random plane wave description, the correlations amongst the Si are entirely negligible. 
This is illustrated in Fig. HI where the average for the covariance is performed over 400 levels 
starting from i = 1500. Both cardioid and stadium billiard exihibit correlations which are 
within the expected statistical errors for being consistent with zero. It is therefore expected 
(and actually turns out) that the semiclassical evaluation of Var[ San] scales as 1/AA^, and 
it is possible to interpret the corresponding multiplicative factor as Var[5j]. 

The remaining task is to evaluate semiclassically the locally smoothed quantity San- 
For this purpose, two ingredients are needed, the wavefunction probabilities |\I'j(r)p and the 
density of particles N(r,E). For the eigenfunctions, a semiclassical orbit summation was 



given by Bogomolny 



32|. It is based on the semiclassical approximation of the retarded 



Green's function (given here for two dimensional systems) 

(31) 



G'^{y,v',E)^^-L= V , ^ exp 



1 TT 

-^^(r,r', E)-i-rif, 



where the sum runs over all closed (i.e. not necessarily periodic) orbits, with the classical 
action of the fi^^ orbit, mi2,^ = dr'j_/dp± the stability matrix element, and t]^ the appropriate 
geometric index (primed and unprimed variables correspond respectively to the initial and 

17 



(a) 1.0 

CoT[SiSi+j] 
0.5 



0.0 



• . . • • 



-0.5 



10 



(b) 1.0 



15 j 20 




FIG. 4: Correlation function Cor [SiSi+j] = Gov /Var [Si] for the cardioid and stadium 

billiard examples. In both cases the correlation function (or covariance) of the billiard is consistent 
with zero to within sample size fluctuations as assumed in the random plane wave model, and as 
implied by the semiclassical theory ahead leading to Eq. (|5ip . 



final coordinates). G^{r,r,E) is related to the local density of states, and thus to the 
eigenfunction probability density via 



u{r, E) = y^ |^^(r)p5(^ - E^) = --Im G^(r, r, E) 

i 

Introducing the density of states p{E) = J2j ~ ^j)^ "we obtain 



(32) 



)a£ 



AN 



(33) 



where on the r.h.s. the overline notation has the same meaning as previously introduced 
except that the division is by AE instead of AA^. 

In the semiclassical evaluation of z/(r, i?), it is typical to distinguish between the "zero- 
length" orbit contribution z/vi/(r) = 171/271 h"^ and the contribution of the remaining orbits, 
whose lengths remain finite as r ^ r'. Here however, interest is in the fluctuations of the 
eigenfunctions, and thus of z/(r, E), near the boundary of the billiard. The orbit responsible 
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for the Friedel oscillations, namely the one returning to its initial location immediately 
after bouncing off the boundary will therefore be extremely short, implying that: i) the 
corresponding contribution will not be sensitive to local energy averaging ; and ii) if r is at 
a distance x from the boundary shorter than, or of the order of, the Fermi wavelength, the 
semiclassical approximation Eq. fl31l) cannot be applied. On the other hand, assuming x 
much smaller than the curvature of the boundary, this contribution can be approximated by 
the exact result valid (in two dimension) for a straight boundary z/Friedei(a^) = ^i^wJoC^kFx). 
This gives 

1 



'^i^)AE = i^iy [1 ± Jo{2kpx)] - -Im G'osc(r, r)^^ , (34) 

valid near the billiard boundary. The tilde on Gqsc indicates that the short orbits giving rise 
to Friedel oscillations have been excluded from the semiclassical sum Eq. (1311) . 

Similarly, the density of states is split into a smooth component and an oscillatory one. 
Once the local energy averaging is performed over a range much larger than the mean level 
spacing, the oscillatory components are small compared with the smooth term. Thus, the 
density of states can be expanded in the denominator. Addressing still two-dimensional 
billiard systems for which pwiE) = Avw gives 

1 1t^^^— 7^ l±Jo(2A;,x) 



^|^,(r)|2^^ = 1 ± Jo(2A;^x) Im Gosc(r, r, E)^^ 7 Posc{E)^j, . (35) 

This equation could be thought of as a slight generalization of the result given by Bogo- 
molny 32], with the only difference that the Bessel function jQ{2kpx) has been introduced 
to account for the Friedel oscillations (which turn out to be important here); see Appendix [Cl 
for an improved normalization of this equation. 

For AE large on the scale of the mean level spacing, but small on the classical scale, the 
energy smoothing can be performed for each orbit contribution noting that dS^/dE = r^, 
with Tfj, the time of travel of the orbit, giving 



with [sinc(a;) == sin(x)/a;]. Energy smoothing therefore implies that orbits with periods 
greater than h/ AE are cut off in the semiclassical sums Eqs. fl3m35p . 

As a direct (and expected) consequence, if the smoothing takes place on a energy range 



larger than the Thouless energy, no orbit can contribute to G'osc(r, r, E)^^^ or posc(-E')^^, and 
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the average wavefunction probability reduces to ^(|\l/i(r)p) = 1 ± jQ{2kpx). Inserting this 
equahty into Eq. f|T6l) with the definition Eq. f|T5l) . we readily obtain the result Eq. f|T7|) for 
the mean value (Si), but here without any assumption regarding the nature of the dynamics; 
i.e. it applies equally well for integrable, mixed or chaotic systems. 



A. Chaotic quantized billiards 



Up to this point, the nature of the dynamics has played no role in the semiclassical 
approach. However, beginning here, the approach is specialized to chaotic systems. The 
oscillating component (a sum over periodic orbits 27|) of the density of states is given by 



Po 



1 1/2 



COS 



h 



(37) 



- 7=pcnodic orbit l^et (M^ - 1)| 

where is the period of the periodic orbit, the monodromy matrix and t]^ the appro- 
priate geometric index. 



For a two-degree-of-freedom billiard, ^|\l'i(r)|'^^^ ^ 1 ± Jo{2kpx) + ef\r)^j^ + e-^''(r)^^ 
with 

2Vh^ 



(1)/ \ 



-Im- 



m 



E 



/^=closed orbit -y/ | i/^i^^l2,/i | 







TT 


exp 


[ h 


^2 



(38) 



AN 



2h 

mA 



(1 ± M2kpx)) 



pert^corbit|Det(M,-l)| 



AE 



(39) 



Note that the /i-orbit sum here includes all returning orbits, not just those in the neighbor- 
hood of a complete periodic orbit. 

One might also be tempted to use the same expressions integrated over energy to deduce 
a similar expression for A^(r; E^). However, the energy integral generates a factor ~ h/TiE 
for the oscillating contribution of an orbit of period r^, and therefore the oscillating terms 
obtained in this way would be of lower order in h than those generated by |\l'j(r)p. In 
leading order, it is therefore only necessary to keep the terms of A^(r; Ef) associated with 
the Friedel oscillations, i.e. begin with Eq. (fT6l) directly and drop further sub-leading terms. 

where 



Thus, Si = S ±. 



c(l) , c{2) 

ijOsc i,osc 



cH _ ^ / Ji(2fc^x) (^) 



(40) 
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(a = 1, 2). Here, two remarks are in order. First, because only the very short orbit contribu- 
tion (i.e the term proportional to Ji{2kpx)/{kpx)) is kept for N{r,E), it is not sensitive to 
the local energy average and can be taken out of the bracket. Secondly, note that the main 
contribution to the integral over space in the r.h.s. of Eq. fHU]) is restricted to the vicinity 
of the boundary. As before we can unambiguously use a system of coordinates r = (x, y) 
with X perpendicular and y parallel to the boundary. To compute S^^^^^, insert Eq. (!39|) 
into Eq. fHOl) . A stationary phase condition has to be imposed in the y direction, but not 
in the x direction since the effective range of interaction is not large, even on the scale of 
the Fermi wavelength. As a consequence, the dominant contributions of the integration in- 
volved in Eq. fl40l) come from the neighborhood of trajectories such that Py = Py (where the 
primed (unprimed) correspond to initial (final) momentum), but for which the initial and 
final x-momenta may differ. Energy conservation however imposes p'^ = ±Px- As illustrated 
in Fig. [5l these trajectories can be associated to a cluster of four orbits which, as x — 0, 
converges smoothly toward the same nearly periodic orbit (or fixed point of the boundary 
Poincare map): i) two nearly periodic orbits such that r lies on the trajectory just before 
or just after bouncing off the boundary, and ii) two non-periodic ones ( p'^ = —px) either 
touching the boundary twice or not at all near r. 

Denoting Si^^i^^y) the action of the nearly periodic orbit to which all of these orbits 
converge as x — > 0, leads to Si{x,y) = ^^^(O,?/) + 6Si{x,y), where 6Si ^ (p^ — Px)x, 
which vanishes for the two periodic orbits and gives ±2|p|xcos6';, with 6i the angle of in- 
cidence of the periodic orbit on the boundary, for the two non-periodic ones. Noting that 
exp(i2/cx cos 6'j) + exp{—i2kx cos 9i) ±2 = 4cs^(A;x cos6'i)^ gives 

sVh . fTiAE\ p , Ji(2A;x) 



*,oscAiv ^ \ 2h J Jn kx ^ ' 

.5K(0,y),(0,y);E) . vr 



(=fixed point 

I f^/^ , 1 

xlm^^= / ay — exp 

V2'ITiJ-C/2 ^y\xl\\x[\\ml2,l\ 



I IVi 

h 2 



(41) 



where the sum runs over all the fixed points of the boundary Poincare section. As in the 
previous section dx '^^ ^^^^^ cs^(kx cos 9i) = [1 ± | sin ^;|]/2/i;. Furthermore, the integral in 
the parallel direction can be performed in a very similar way as in the derivation of the 
Gutzwiller trace formula. Using the fact that near the periodic point (0, yi) 

S,m,m + Sy), {0,m + Sy); E) = S,{E) + D^t(M,-l) ^^, ^^2) 
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(a) 




(b) 




(0 




(d) 



FIG. 5: Sketch of the four orbits which, as x — > 0, coalesce into the same [nearly] periodic orbit. 
The top row corresponds to two (nearly) periodic orbits such that r lies on the trajectory either 
(a) just before or (b) just after the bounce off the boundary. The bottom row corresponds to two 
non-periodic orbits (p^ = —px) such near r (c) one of them does not touch the boundary and (d) 
the other one touches twice. 
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f F I COS 6i I we get 
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SijE) 
h 



r, 



j,osc^jV 



:Sinc 



2h 



(43) 



with \{9) = [1 ± I sin6'|] /I cos6'| the same function that was introduced in the random plane 
wave approach (cf. Eq. fl26|) ). 



The computation of 5^ olcAiv sven simpler as the only spatial dependence of ef'^ arises 
from the Bessel function Jo(2A;^x). To facilitate the comparison with Eq. ( H3l) . replace the 
sum over periodic orbits by a sum over fixed points of the Poincare section, in which case 
the period of the periodic orbit has to be replaced by the average time of flight r^/n^ = 
i^/vpn^ {i^ and are respectively the total length and total number of bounces of the 
peridic orbit 7) between two successive bounces on the boundary. Using J^{du/u) Ji{u){l± 
Jo{u)) = (1 ± 2/7r) gives 



(2) 

ijOsc^jV 



a. 



Akl ^ IniA 

t=nxea point 



1 ± 



COS 



h 



r, 



vr; ^1 Det(M, - 1) 



smc 



TiAE 

2h 



(44) 
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def 

where ni = n^{i) is the number of bounces of the periodic orbit to which I belongs, and 
a similar notation is implied for the other parameters. This actually depends only on the 
periodic orbit 7 and not on the specific periodic point I on 7. The (2/7r) factor can be traced 
back to the jQ{2kpx) term, and therefore eventually to the Friedel oscillations of the local 
density of states. 

The two terms can be combined to give 



— u. 



7 2 



'Akl 



7=periodic orbit 
Z=fixed point of 7 



v/| Det(M 



smc 



m 



7 



T^AE 

2h 



2Ani 



71 



(45) 



This form is suitable for performing the calculation of the variance. However, as rederived 
briefly in Appendix [Bl note that the mean length per bounce in a billiard is 



d 



(46) 



Also, the density of fixed points is uniform for long orbits in the measure d sin 6. Averaging 
with that measure while replacing Cii/2Ani by its mean value 7c/2 gives 

'.7r/2 



2Ani 



1 ± 



7r 



dsin^ 



-tt/2 



(1± 


sin 


Oil) 




COS 6i 





2 



1 ± 



TT 



2 2 



1 ± 



TT 



= (47) 

Thus, the constant (l ± ^) can be understood as arising from ^^^osc^^ (i.e. associated 
with the density of states) as the angular mean (A(6'))e of the corresponding term in S^^^^^^^. 



To compute the variance, it is necessary to square Eq. fl45|) and average the re- 
sulting expression over a large energy range. For two periodic orbits 7 and 7', 



cos 



v - 

2 



COS 



h 



r, 



equals one half if 7 and 7' are either the same or- 
bit or time reversal symmetric, but zero otherwise, which makes cross-terms from different 
periodic orbits vanish. Note however that it does not eliminate cross-terms of the various 



fixed points for a given orbit since these contributions oscillate with the same frequency. 
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This gives 
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(48) 



For long orbits, which are going to dominate this sum, it is possible to identify i.y/n^ with 
d = ttA/C, the mean length per bounce in the billiard, and assume that the angles 61 are 
uncorrelated and uniformly distributed with the measure dsin6'. It turns out for the last 
sum in Eq. ( HHl) 
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(49) 



where Eq. fHTj) has been used to cancel the cross terms between different fixed points. Making 
use of the Hannay-Ozorio de Almeida sum rule 41| in the form 

1 



E 



fixed points I 
with n bounces 



Det(M, 



1, 



(50) 



(where the sum runs over all fixed points belonging to a periodic orbit with n bounces), 
identifying the period r of the orbit with nd/vp {vp = hkp/m is the Fermi velocity), replacing 
the sum over the number of bounces by an integral, and making use of Eq. ( l46l) gives 



{{s)AN-sy 



2 k^A^ 
9 s kpC 



{{\{e)-{X),y)e / dnsinc 



\X{e) - {X)ef), 



ndAE\ 
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(51) 



2 27r3 \^ ^ ' N i-j le 
As expected, the variance of S^n is inversely proportional to AA^ = pw{E)AE. Applying 
Eq. (l30l) . the absence of a term constant with AA^ confirms that, as assumed in the random 
plane wave approach, there are at this level of approximation no correlations amongst the 
Si. Thus, the variance of the Si are given by 

kp/l: 



Yar{Si 



27r3 



X ((A(^^) - {X),r)e 



(52) 
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f) 



Dirichlet 



Neumann 



(53) 



B. The periodic orbit spectrum 

Beyond the Var[iSj], which here characterizes the universal (long time) behavior of the 
system under consideration, the semiclassical treatment developed in the previous subsection 
provides information on system specific quantities. In particular, it makes it possible to 
address phenomenon related to shorter time dynamics, and thus quantum mechanically, to 
longer energy range. For example, Eq. (H5l) can be used directly to compute the Fourier 
transform of the Si. Interestingly, Eq. ( H5l) has a structure very similar to that of the 
density of states Eq. fl37|) . This gives a simple and striking prediction, namely that the 
Fourier transform of the Si will display peaks at the same locations and with the same 
shapes as the Fourier transform of p{E) (up to the transformation i ^ E), and will be 
simply scaled by factors which depend only on the lengths of the orbits and on their angles 
of incidences {61} at the various places where they bounce along the boundary. In Fig. [6](a), 
the Fourier transform of the cardioid billiard density of states is shown in comparison with 
the Fourier transform of the {Si} displayed in Fig. [6](b). The peaks are in precisely the 
same positions and their shapes are similar, but the amplitudes differ as expected; as a 
parenthetical remark, for technical reasons the Fourier transform of the {Si} uses a slightly 
different Fourier transform than the density of states (effectively divided by the wave vector), 
which is denoted by a subscript in the remaining figures, but this has no effect on the overall 
discussion of the physics involved. In addition, the prediction of Eq. ( H5i) for the shortest 
periodic orbits and their retracings is shown. The predicted amplitudes for the {Si} are 
reasonably close, although perhaps slightly too large by 30-50%. Otherwise, the agreement 
with Eq. fHSjl is excellent. The excess in the prediction for short orbits is curious because 
if the predicted amplitudes for all of the orbits were too large, it should be found that the 
prediction of the variance is slightly too large instead of a bit too small (say factor of two) as 
in this case. We have checked that in fact, the predictions for long orbits, which dominate 
the calculation of the variance of the {Si}, are indeed a bit too small, the opposite of the 
short orbits. Why it has turned out this way for this particular example remains for future 
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FIG. 6: Fourier transform of the density of states and Si — {Si) of the cardioid billiard using 
Dirichlet boundary conditions. In (a), the solid line is for the density of states using the first 2000 
odd parity levels of the cardioid billiard. The dashed line is for the density of states using the 
corresponding form of the Gutzwiller trace formula. In (b), the solid line is the Fourier transform 
of Si — {Si) for the first 2000 Si and the dashed line is the result of Eq. (j45p calculated for the 
shortest periodic orbits as the one shown in Fig. [T] and the insets (along with their retracings). 
The agreement is quite good. 

consideration. 

As a last remark, note that the tendency for long orbits to explore uniformly the phase 
space implies both that sin^; is distributed uniformly and that the mean length between 
bounces l^jn^y for a given orbit 7 can be identified with rf, the full system average distance 
between bounce. This is what made it possible to apply Eq. (H7|) and to cancel the cross- 
terms between various fixed points of the same orbit in Eq. fH^ . Had the term proportional 
to in Eq. fl49l) not been zero, it would have given rise to a contribution parametrically 
larger (in K) than the computed one [~ {kpCY instead of {kpC)]. Short orbits, for which the 
cancelation of cross-terms done in Eq. fHOl) cannot be applied may therefore have a stronger 
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influence on the fluctuations of tlie iSj's tlian wliat miglit be naively expected from Eq. (15T!) . 



C. The basic distinctions in the two theoretical approaches 

Interestingly enough, the expressions in Eq. fl53|l are exactly the same results as the 
random plane wave approach Eqs. (12711281) . except for two differences. First, the mean square 
{X'^{6))g has been replaced by the variance of X{0), giving now a much better agreement with 
the billiard results (see FigO]). Second there is a factor two difference in the prefactor. 



1. Proper normalization 

The replacement of the mean square by the variance can be related to the lack of proper 
normalization of the wave function in the random plane wave model. Indeed, fluctuations 
of the local density of states z/(r) can either imply fluctuations of the wavefunction proba- 
bilities |\E'(r)p, which have to integrate to zero because of the wavefunction normalization, 
or fluctuations of the total density of states for the part which survives the integration over 



space. The role of the term proportional to posc{E)^^ in the right hand side of Eq. ( l35i) can 
therefore be understood as ensuring the normalization of the eigenfunctions. As this term 
is precisely the one giving rise to the contribution proportional to (A)^, it turns out that 
in the semiclassical calculation, proper normalization of the eigenfunction is what generates 
the variance of X{6) rather than its mean square. Since the random plane wave model used 
here imposes normalization on average rather than for each individual eigenfunction, this 
contribution is necessarily missing there. A modifled version of the random plane wave 



model in which normalization is better enforced 
this issue. 



43( 1 should, however, properly address 



2. Dynamical correlations 

The factor two difference in the prefactor (or conversely, the fact that except for this factor 
two and the normalization effect, the random plane wave and the semiclassical expressions 
are identical), although less important from a quantitative point of view, is however puzzling 
enough to deserve further discussion. To focus better on the main point, consider two 
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simplifications of the problem under discussion. First, assume as understood the issue 
of eigenfunction normalization, and consider below only the contribution from the Green 
function (i.e., ignore density of states fluctuations). Second, consider that the procedure 
used to extract the variance Var[iS.j] from the locally smoothed quantity S^n [see Eq. fl5U]) ] is 
equivalent to the effective rule according to which the various quantities under consideration 
should be smoothed over an energy window of width A (so that AA^ = !)• 

Having this local smoothing in mind, and ignoring for the moment the fluctuations of the 
density of states (i.e. assuming there is exactly one state in each interval 6) gives 
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Close to some reference point r and not considering yet the proximity of a boundary, this 
gives for the oscillating part of the wavefunction probability 
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Above, the sum runs over all closed trajectories /i starting and ending on the reference point 



r, with initial and final momenta p^ and p^, time of travel r^, and 



A, 



1 



1 



ih 2mh 



exp 



' h 



TT 

1^2 



(57) 



In the semiclassical calculations of section [TV] it is taken into account that as the integra- 
tion over space is performed, closed trajectories are continuously deformed, and in particular 
the initial and final momenta p^ and p^ vary. As a consequence the dominant contributions, 
which correspond to nearly periodic trajectories (to within a bounce off the billiard bound- 
ary in this particular calculation), can be understood as arising from the neighborhood of 
periodic orbits, leading to the periodic orbit sum Eq. fj43|) . The calculation of the variance 
is then done using the Hannay-Ozorio de Almeida sum rule Eq. (ISUjl . 

Consider that if the dynamical correlations, i.e. variations of the orbital properties (initial, 
final momenta p^, p^, prefactor A^, etc..) are neglected, the semiclassical expression for 
vl>*(r')\l>(r') greatly resembles the random plane wave model. Indeed, if long orbits are 
dominant: 
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i) initial and final momenta and are independent and uniformly cover the energy 
surface; i.e. the model can be taken as a "random pair of plane waves" model. 

ii) applying the diagonal approximation in the semiclassical calculation amounts to 



A*^A^i oc 6f^f^i or, if the system is time reversal invariant /i and fi' are related through time 
reversal invariance. Near a boundary, the correlations are included between the trajectories 
related to one another by a bounce off the boundary. 

iii) although the are not Gaussian distributed, the fact that the number of trajectories 
is extremely large for long orbit makes it possible to use a central limit theorem, implying 
that only the variance of these quantities are relevant (and that one can as well consider 
them as Gaussian). 

iv) the variance of the is constrained by the sum rule valid for closed orbits (a slightly 
different rule than the Hannay-Ozorio de Almeida sum rule used for periodic orbits) 4J], 

^ \A^\^S{t - r^) = —uwPciir, r, r) (58) 

where for long orbits in billiards, the probalitity of return Pci(i", r, t) can be taken uniform 
and equal to 1/^. Properly carrying out the smoothing on the range A produce the damping 
factor sine (^^f^ j of Eq. (l36l) . so that 

Thus, neglecting the spatial variations of the orbits properties which contain dynamical cor- 
relations, quite standard semiclassical approximations, namely the diagonal approximation 
and the assumption that P(r, r, r) is uniform, makes it possible to derive a "pair of random 
plane waves" model, not completely identical to the original random plane wave model, but 
similar in spirit. It can be shown furthermore than computing Var[5'j] under this model 
gives exactly the same result as the random plane wave model. 

Indeed, one can compute (ei(ri)ej(r2)) - (ei(ri)) (ei(r2)) = A'^ {\iJi{ri)\l^^\tpi{r2)\lsc) in 
this way [see Eq. (12(1 )]. To be more precise, let the {A^} be uncorrelated unless the 
corresponding trajectories are time reversal symmetric, or related one to each other by a 
bounce off the boundary of the billiard near the initial or final point of the trajectory (i.e. 

((pD/.'>(p^)m') = (±(pi)M' (Py)^) and {{pD , (pD ^,') = (±(p^)^, (p^) J . If the reference 
point r is taken on the boundary, and measuring the distance x from the boundary, this 
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amounts to taking equal for these trajectories, giving 
((l^^(ri)P).JI^.(r2)r)„J = 

C^^Wy - kl)iyi - y2)]cs{plxi/h)cs{plXi/h)cs{plx2/h)cs{plX2/h) , 

(60) 

which, using the sum rule Eq. fl59p and inserting the resulting wavefunction correlations in 
Eq. flTS]) gives exactly Eq. fl22]) derived with the random plane wave model. 

To summarize, neglecting dynamical correlations in the semiclassical approach generates a 
random pair of plane waves model that, in essence is derived with usual approximations. For 
the problem we consider here, this model gives exactly the same result as the random plane 
wave model. The "random pair of plane wave" is however not ad hoc, whereas the random 
plane wave model is. This makes it possible to discuss precisely what approximations have 
been made, and therefore in what way we could expect the random model to differ from the 
purely semiclassical treatment. In particular we see that the interferences between reflected 
wave at the boundary is treated in the same way in both the semicalssical and the random 
approaches, giving rise to the same \{9) dependence. On the other hand the prefactor is 
related, in the semiclassical approach, to the way classical orbit with nearly matching initial 
and final momenta are structured around periodic orbit. This aspect is completely ignored 
in the random models, and we therefore cannot expect them to give exactly the correct 
prefactor. 



V. NON-CHAOTIC SYSTEMS 



Included in the class of non-chaotic dynamical systems are three main subclasses: i) the 
limiting case of integrable systems, all of whose dynamics are regular; ii) near-integrable sys- 
tems, characterized by having classical perturbation theory generally work well in describing 
its dynamics; and iii) mixed systems, which contain an intricate mixture of both regular and 
chaotic dynamical regions in their phase spaces. Generally speaking, semiclassical theories 
and what is known vary according to each subclass. For example, trace formulae exist for 
integrable 



46l | and near-integrable systems [47|, |48|, but not for mixed systems. In fact. 
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a proper treatment of semiclassical theory for mixed systems is lacking. However, for the 
purpose here of investigating the properties of the set of {Si}, only the simplest level of 
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semiclassical theory is considered. In other words, in regular dynamical regions (whether 
from integrable, near-integrable, or mixed systems), structures called tori are assumed to 
exist, which are invariant manifolds under classical motion, and possible complications from 
resonances, diffraction, or tunneling are ignored. In chaotic regions, only the complication 
of a family of marginally stable orbits is considered beyond that which was already treated 
in the previous section. 

Unlike chaotic dynamical regions in phase space, for regular regions there are two possi- 
ble overarching semiclassical approaches. In the first, particular tori quantize allowing the 
detailed evaluation of Sj for each eigenstate. This is based on the Einstein-Brillouin-Keller 
(EBK) scheme 49|, |50|, l5l|]- In the second, a periodic orbit trace formula results from ap- 



plying a Green function approach, much like the chaotic case. However, this Green function 
approach is not given here since the information about each individual Si, the mean, and 
variance are already understandable through the EBK approach. 

A. Einstein-Brillouin-Keller quantization 

1. General expression for Si 



Continuing with two dimensional billiards, each torus is characterized by two action 
variables (Ji, J2), and it is always possible to choose the corresponding angles {ipi, (p2) such 
that the intersections of the torus with the boundary of the billiard are parameterized as 
V^i = /k(v^2)! = 1, • • • , /^max, with Kmax the uumbcr of bounces on the boundary for the 
considered torus. 

An eigenstate is constructed on a quantizing torus 



Ji = 27r7i(n[*-* + cri/4), J2 = 27r/i(n2'' + 0-2/4) where (cri,cr2) are the Maslov indices 



and can be expressed as 



27r 



E 



5(V5l,V22) 



d{x,y) 



exp ( jSe,{x, y) 



(61) 



The sum runs over the various sheets of the torus projecting onto the point r = (x, y) and 
Si{x,y) is the corresponding action (including the Maslov phases). 

Consider in greater detail, the neighborhood of the ipi = /k(v52) boundary. To further 
simplify the discussion, assume that the torus has only two sheets (corresponding to negative 



31 



and positive ipi — ft^{ip2)) projecting on any given point (x, y) near this boundary. The results 
derived under this hypothesis apply in the general case, as is justified below. Adding the two 
Vi < /k(v^2) and ipi > /k(v^2) contributions, and expanding the action from the boundary 

as S{x,y) = S{x = 0,y) +PxX, generates 



exp ( -S{x = 0,y) ] 2cs{pxx/h) ; 



(62) 



d{x,y) 

assuming that the local variation of the Jacobian determinant in the direction perpendicular 
to the boundary can be neglected. Inserting this expression into Eq. f|T6l) we obtain 

C 



with 



5,; 



S,. 



dxdy 



It 



— cs {pxX/h) . 



(63) 



(64) 



d{x,y) 

If the torus has more than two sheets projecting onto the neighborhood of the boundary (in 
which case Eq. (!62|) involves a sum), the rapidly oscillating phases, {exp{iS{x = 0,y)/h)}, 
eliminate cross terms upon integration over y, and thus the calculation of Si would involve 
just a single sum. 

Changing the integration variables to (y)^, t), with r measuring the time from the bounce 
on the boundary of the billiard and the angle ip2 at that bounce, we can further simplify 
Eq. Indeed 



y^2 = f2 + 

= /«(<^2) +^1^ 



(65) 
(66) 



with Ui = dH/difi (i=l,2) the angular frequencies, and therefore the Jacobian can be 



expressed as J 



\uJi — LJ2{dfK/dip2)\. Noting furthermore that x{{p2,T) 



Vp cos{6{ip2))T, the integral on the variable r can be performed explicitly, giving 

/^max -I /*27r 

Si = L - u;2{dfjdy,'2)\\ [OM)] ■ 



K=l 



(67) 



2. Circular billiard 



Beyond the orders of magnitude, the explicit evaluation of the variance of the set {Si} 
for integrable systems is very much system dependent as it is usually not possible to 
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make general assumptions about the correlations between the various quantities involved 

(A( Ji, J2, 6'k(v22))? '^i,2(</i, </2), c^/K/rfv^^l"^!' "^2, 6'k(v^2))- We therefore consider now a specific 
system, namely the circular billiard. 

The computation of the expression Eq. fpTl) for this billiard is made somewhat simpler 
because the angle 9 at which trajectories bounce off the boundary is a constant for a given 
torus, and thus a function of the actions (Ji, J2) only. Furthermore, for a given invariant 
torus, we can construct the paths on which the action variables ( Ji, J2) are constructed as 
the cut of the torus in the radial direction (for Ji), and the caustic (or any topologically 
equivalent path) for J2. In this way, the angle (p2 can be identified with the angle of the 
polar coordinates, and the boundary of the billiard can be taken as <y9i = (i.e. = 
0). Furthermore the angular frequency ti;i(Ji, J2) can be identified with 27r/t(Ji, J2), with 
i{,Jii J2) the time between two successive bounces for trajectories of the corresponding torus. 
The expression Eq. (!67|) thus takes the simple form 

^ 2i A[g(Ji,J2)] _ .gg. 

kpVp t{Ji,J2) 

Quantizing tori sample uniformly the plane (Ji, J2), and therefore statistical quantities 
such as average and variance should be computed with the measure dJidJ2. However, 
using that the change of variable (Ji, J2, v^i, ¥^2) {E,^,T,p^) introduced in appendix [Bl 
[p^ = Pp sin 6) is canonical, implying dJidJ2dipid(p2 = dEdp^drdC,, and that for the circular 
billiard 6 and E depend only on the action (Ji, J2), and r and ^ on the angles {(fi, (f2), one 
can write, in the neighborhood of the Fermi energy Ep 

dJMSiE - Ep) oc t(Ji, J2)d{sm9) , (69) 

and thus use the probability measure 

with {T{6)) given by Eq. ( 1B6I) arising from the normalization (given for allowed values of 
sin 6* in the range < sin^^ < 1). 

One can in this way recover in the particular case of the circular billiard the general 
expression Eq. f|T7|) of the mean value (Si). Indeed 

2? r 

(5,) = ^^(A(^)). and 

{Si) = iilT 7^ I ± "^i^—r (^±1) = i(l + -^^] . (71) 
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Similarly, the expression for the variance reduces to 



Var [Si 



k^VpTlA 

klnR^ 



'd(sin0)^- 
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(72) 



cos Q cos^ ^ 

As before, the constant for Dirichlet boundary conditions is quite small, and the divergence 
of the Neumann boundary conditions increases the order. In this case, the effect is greater 
than logarithmic. Using the same cutoff for the Neumann case as at the end of Sect. IIIIB 
and that i ~ Ak^j^-n^ we obtain 

, 2 



Var [Si 



I X 



^ - i - (1 - ^ 0.00457 Dirichlet 

7r 2 V TT / 

_2 _ 1 _^ ^1/2^1/4 _ (1 _ Neumann 



(73) 



Thus the variance scales proportionally to i for Dirichlet boundary conditions and t'l'^ for 
Neumann boundary conditions. 



B. Bouncing ball modes in the stadium billiard 

Even for fully chaotic systems, it is possible to have a situation where some (with vanishing 
measure) of the trajectories behave more like those of integrable systems. An example is 



provided by the bouncing ball orbits of the stadium billiard [52(]. Tanner [53|] showed that for 
the purposes of a semiclassical theory of eigenstates, the phase space in the neighborhood of 
the bouncing ball orbits behaved much like an island of regular motion, and that families of 
orbits that cannot be taken as isolated contribute in essential ways. This greatly complicates 
the desire for a rigorous semiclassical theoretical approach. Though these states can be 
thought of as EBK-like states similar to those studied in Sect. IV Al they do get connected 
through diffractive terms to the chaotic states, and thus some of them behave more like 
resonances rather than individual quantized states. 

A complete semiclassical description is beyond the scope of this study. Indeed, as men- 
tioned in the beginning of Sect. |Vl as many dynamical system complications as possible are 
being neglected here. Instead of attempting a rigorous semiclassical theory for the bounc- 
ing ball modes, a rough approximation is given instead. To be specific, we consider here 
the even-even symmetry states of a stadium billiard with Dirichlet boundary conditions, or 
equivalently a symmetry-reduced quarter stadium with Dirichlet boundary conditions on 
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the original boundaries and Neumann boundary conditions on the symmetry hues. To start, 
consider the eigenstates of a rectangle with one side length equivalent to the side length Lg 
of the symmetry-reduced quarter stadium and the other, the radius of curvature R. These 
states can be used to give an approximation to the bouncing ball modes. In essence, the 
bouncing ball modes with few nodes along the side length (ignoring mixing into the chaotic 
states) vanish quickly upon entering the quarter circular end-cap. A quantization along the 
side length direction with Dirichlet boundary conditions on the side entering the end-cap 
and along the side length itself is a good starting point. Since our calculations have been 
done for even-even symmetry eigenstates, consider Neumann boundary conditions for the 
remaining two sides. The normalized states are given in cartesian coordinates by 



\^i{qi,q2)\' 



■ cos 



2m + 1 



Tcqi cos 



2n + l 
2R 



-vrg2 



(74) 



where the origin is the corner, m is a small integer, say 0,1,2,3 or so, and most of the 
kinetic energy is in the other direction and so n is a large integer. 

To write the equation for the Si, the coordinate system of the state must be rotated and 
translated to the boundary coordinate system used for the Friedel oscillations separately 
along the two symmetry lines and the side length. This gives three terms to evaluate for the 
bouncing ball contributions, 

Cd-Cn\ . 2i f°° ^ Jii2kpx) ,/2m + l 
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Cd~Cn + 2A\- + — 



(75) 



TT \ R Lg 

where the straightforward integrals over the y coordinates have been evaluated before writing 
the first expression, and Cr, = + tiR/2 and = Lg + 2R are the perimeter lengths with 
Dirichlet and Neumann boundary conditions respectively. In the last expression, the overall 
mean Eq. ( |T7I) is subtracted and the square root reduces to unity since bouncing balls with 
significant momentum toward and away from the end-cap do not exist. 

There are a number of interesting consequences of Eq. (!75l) . The last line captures the 
scale of the deviation from the mean. Putting in the parameters used for the stadium 
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calculations of Figs. (12][3!) {R = Lg = I) generates a growing expected deviation from the 
mean that hits about 30 for i = 2000. This result also implies that Si for the non-bouncing 
ball modes fluctuate about a negative bias given by 



^(non bb) 



bb 



/bb kpA 
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7T 



2A 



1 

R 



1 



(76) 



where the fraction of bouncing ball modes is denoted /bb- The best known estimate of 
/bb for the stadium billiard to our knowledge comes from Tanner jsS], which implies that 

3/4 

i we have introduced the ratio 7 = Lg/R and used that Am = Akp. 



/bb ~ 7 |^^(^+^/4) 

For other systems with bouncing ball modes, their fraction may scale differently [5^. The 
results of Eqs. f l75H76l) are given by the dashed and solid lines in Fig. [2](b). The dashed 
(upper) line crosses right through the neighborhood of the peak values. For such a rough 
approximation, Eq. (1751) is quite good. In addition, the solid (lower) line captures the 
negative bias implied for the non-bouncing ball modes quite well. 

Second, to the level of approximation here, any of the bouncing ball modes contribute 
fairly equally locally in energy (determined by the quantum number n); i.e. there is almost 
no dependence on sequence number m. Although, the peak values from the bouncing balls 
do not appear to be constant in the figure, presumably due to weak admixtures of chaotic 
states, this feature greatly simplifies a calculation of the fiuctuations. The variance can be 
inferred by noting the following: i) the square of the bouncing ball deviation from the mean 
multiplied by their fraction denoted /bb gives their contribution; and ii) the contribution of 
the remaining 1 — /bb non-bouncing ball modes can be taken as the square of the average 
amount they must each be in deficit of the mean plus the [sub-leading] variance from the 
chaotic system results. The combined consequences lead to the expression 
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(77) 



Taking account of the decreasing fraction of bouncing ball modes with increasing i, the 
variance scales as i^^^. See Fig. [3] for a comparison of this formula with Var[5j] for the even- 
even eigenstates of the stadium billiard. Again, the rudimentary approach here captures the 
main behavior fairly well. 
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FIG. 7: Fourier transform of Si — {Si) and the density of states from the even-even stadium 
eigenstates using Dirichlet boundary conditions. The dashed hne is for the density of states (which 
has a different vertical scale shown at the right side). 

To conclude this subsection on the bouncing ball orbits, a few remarks are in order. 
First, we stress that although the classical dynamics of the stadium billiard is, mathemat- 
ically speaking, purely chaotic, as far as Si statistics are concerned the existence of the 
marginally stable bouncing ball family makes this system behave very much like an inte- 
grable billiard. In particular the scale of the fluctuations are order of magnitude larger 
than for the "genuine" chaotic billiard considered in sections [Till and HVl Second, because 
of the presence of two classes of states with drastically different properties, the covariance 
amongst the Si is, again in contrast with genuine chaotic systems, non zero. Furthermore, 
both the variance and the covariance show an energy dependent structure which complicates 
significantly the extraction of these quantities from the locally smoothed iSjaat as was done 
in section HVl 

As a final remark, although the periodic orbit formula is not derived here. Fig. [7] is shown 
for completeness. As in Fig. El the results are compared between the density of states and 
the {Si}. Again, the peaks are in the same positions, i.e. those determined by periodic 
orbits, but with differing amplitudes. The importance of the bouncing ball modes is quite 
visible. 



VI. DISCUSSION 



The mean and variance of the quantities {Si} treated in 25|] and in greater detail in 
this paper have been introduced as examples of a new, and non-local, class of statisti- 
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cal measures with physical relevance; the {Si} are connected to the addition spectrum of 
quantum dots. By no means should they be taken as the only possible measures repre- 



senting this c 
Fermi gasses 



ass. Indeed, in recent works on finite-size fiuctuation properties in ultracold 
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551], involving fiuctuations in the Bardeen-Cooper-Schrieffer pairing gap. 



a quantity was introduced for the order parameter that is similar to the {Si}, though with 
additional complications. Not surprisingly, several common features are found for both 
quantities. As the fiuctuations are dominated by a term arising from the interplay of the 
Friedel and eigenstate oscillations, there is a significant decrease in fiuctuation magnitude 
due to Dirichlet boundary conditions, either through reduction or even vanishing of the 
prefactor of the leading term in kpL. In the case of the {Si}, the prefactor is decreased 
by more than an order of magnitude with respect to the prefactor for Neumann boundary 
conditions. Another critical feature is the role of dynamics in the scale of the fiuctuations. 
Chaotic dynamical systems lead to a fiuctuation scale of lower order in kpL than integrable 
or mixed dynamical systems. Again, this leads to a fiuctuation scale decreased by an order 
of magnitude or more. The decrease in scale for chaotic systems can be traced to an ergodic 
nature of the individual eigenstates. Conversely, the much larger fiuctuation scale for in- 
tegrable and mixed phase space systems, suggests the possibility of new physics associated 
with more regular dynamics or the possibility of measurements that can be used to deduce 
information about the dynamics. 

The i^/^ oc kpA/C (or i^/^lni) dependence for chaotic systems, and faster-growing de- 
pendence for integrable systems, implies that the fiuctuations embodied in Var[iSj] in fact 
grow with the size of the system, becoming eventually larger than of order unity. This 
implies that the corresponding residual interaction contributions will in that case become 
larger than the mean level spacing A. In other words, since A is the energy scale set by the 
one particle energies, the fiuctuations in the residual interactions may become large enough 
that they generate a modification in the ground state orbital occupation number, and more 
generally reach a point where a first order perturbation treatment of the interactions is not 
adequate. 

Returning to the two theoretical approaches included here - a random plane wave model 
for chaotic systems and semiclassical theory, whatever the dynamics - we have seen that the 
basic random plane wave model is intrinsically less powerful than semiclassical theory, but on 
the other hand, it is technically much simpler to implement. Surprisingly, given the excellent 
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results it generates in other contexts, the random plane wave model here displays a couple 
of significant faults. The most simple to track down is the effect of having its normalization 
be across the ensemble as opposed to the individual eigenstates. For chaotic systems, this 
led directly to the replacement of the mean square of A(^^) by its variance. Inclusion of the 
variance improved the theory considerably as the mean square largely overestimated the 
prefactor. In addition, the proper treatment of dynamical correlations in the semiclassical 
theory led to a factor two increase in the prefactor constant. Truth be told, for the cardioid 
billiard example, the results seem to agree better without the factor two (see Fig. [6]), but a 
concerted search for an error in the semiclassical calculation never resulted in its removal. 

Finally, it is important to be aware of some consequences of decompositions such as given 
in Eq. The separation of average and fluctuating parts of a non-local statistical measure 
may not be the full story from a physical perspective. In fact, this is the case for the {Si} 
treated here. As noted earlier, the mean of Si does not lead to any modifcation of the ground 
state occupations numbers as its effecs gets cancelled. On the other hand, the fluctuating 
part of Si does affect the ground state, but there are two components. The leading order 
fluctuations that come from the use of the term Ai"sec(r; E^) in the fluctuation expressions is 
essentially a mean fleld effect, analogous to scrambling. These fluctuations, if large enough, 
can reorder the fllling of the single particle levels. They do not however lead to high spin 
states or other unusual behaviors. The remaining term 6N{r] Ef) is responsible for exotic 
physics in those cases where it is sufficiently large. As the focus throughout this paper 
was on developing the theory of the non-local statistical measures themselves, the leading 
behaviors have been emphasized, which though dominant, are not necessarily the only ones 
that deserve to be considered - that depends on the physical context (other types of problems 
exist, such as the fluctuations of superconducting gap mentioned ealier, where the dominant 
terms in the relevant non-local statistical measures do contain all the important physics). 
Therefore, revisiting the fluctuation-fluctuation term involving 5N{y] Ef) will sometimes be 
important, but is left for future work. 
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APPENDIX A: SECOND ORDER TERMS IN THE AVERAGE PROPERTIES 
OF THE Si 



There are second order terms in the computation of the average properties of the Si 
coming from the boundary, which are given here. Although, the other second order terms 
from the curvature and boundary discontinuities are not being derived and hence this cal- 
culation is incomplete, the numerical calculations necessary to isolate the average behavior 
before calculating the variance or covariance are improved by including them. Therefore, an 
account is given here. 

The decomposition of Eq. giving 6N{r] E^) implies 

SN{r;Et) = Y.mr)f-{\^,{r)f) . (Al) 



Substituting the relations from Eqs. fl9l[T5|) and integrating the constant terms gives 

S, = i±-!dv —^^-^^^^e.{v)±A5N{v-E)e,{v) . (A2) 

The last term merits some discussion: its leading behavior is seen to be two orders weaker 
than the overall expression and is straightforward to evaluate because only its leading contri- 
bution is required. Under the operation of taking the expectation value, the only surviving 
term is 

(5Ar(r;E)e,(r))=^([|^.(r)|^-(|V^,(r)|^>]') , (A3) 

so that after integration over space, this is essentially equivalent to {Ma) — (Mij). Indeed, 
the local Gaussian random behavior is uncorrelated from state-to-state in the random plane 
wave model so that the only surviving term comes from the i*'' state with itself, all others 
vanishing. Locally, before squaring and taking the expectation value, the expression on 
the r.h.s. of Eq. flA3l) can be thought of as being like the square of a zero- mean, unit- 
variance Gaussian random variable with unity subtracted. However, it does have a variance, 
which is position-dependent and given by the right-hand-side of Eq. f[T^ : i.e. the r.h.s. acts 
as an envelope. Inside of the billiard (excluding the semiclassically vanishing boundary 
region), its value is however constant and equal to the inverse of A to leading order. This 
generates a constant, equal to two, after integration. Therefore, the expectation value of Si 
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is approximately 
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where we choose kp = ^j2mEi/h. Interestingly enough, if a length of the boundary 
follows Dirichlet conditions and a length Cn follows Neumann conditions, it is not correct 
just to make the substitution ±£ ^ Cn — Cd- Rather, the above expression becomes 
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after redoing the algebra. The distinction arises because some of the correction terms depend 
on the sign of the boundary conditions, whereas other correction terms depend on the sign 
squared. 



APPENDIX B: MEAN LENGTH OF A TRAJECTORY BETWEEN TWO SUC- 
CESSIVE BOUNCES 



This appendix briefly rederives Eq. P6|) . which gives the mean length of a trajectory 
between two successive bounces off the boundary of a billiard. While this is a well-known 



result (see 56] and references therein), several equations used in the derivation are needed 
in the main text. The result can actually be obtained by computing in two different ways 
the energy surface volume of the billiard 

V{E) = J dpdrS{E - H{r, p)) , (Bl) 

with H{p,r) = p^/2m + V"(r), and V{r) = inside the billiard and oo outside. The first 
way is to perform this integral with the original coordinates (r, p), giving V{E) = 2iTmA. 

The second way to perform this integral is to use another set of coordinates, constructed 
as follows. Any point (r, p) of the billiard's phase space, can be considered as belonging 
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to a trajectory which has last bounced off the boundary a time r ago at a location on the 
boundary labeled by the curvilinear abscissa ^. Denote ro(^) the corresponding point on the 
boundary and introduce the action 

^(ri,r2,r,0 = / L{r')dT' (B2) 



with L = pr — H the Lagrangian function. Since dS/dr = p, S{ri,r2,T,^) can be used as 
the generating function of the canonical transformation 

(r,p)^(Q,P) (B3) 

with Q = r). The new momentum coordinates are thus given by 

P - dS _ dS _ 

^ _ OS _ dS _ 

with the projection of the momentum p on the direction parallel to the boundary at ^. 
As Eq. ( ]B3I) is a canonical transformation, drdp = dQdP and the energy surface volume 

is 

V{E) = J d^dp^dEdrSiE - H) . (B4) 

Performing the straightforward integration over energy, and noting that for given {C,,p^) the 
integral over dr yields the time of travel t{^,p^) of the corresponding trajectory between ^ 
and the following bounce gives 

V{E) = j d^dp^t{^,p^) = 2pCt , (B5) 



with p = \/2mE and t the mean time of travel between two successive bounces. Identifying 
this expression with the one obtained using the original coordinate system gives 

^* = ^ (B6) 

m L 



and finally, using pjm = |r|, one finds Eq. (j46 



APPENDIX C: NORMALIZATION CORRECTIONS TO EQ. 



Note that in Eq. (139|) only the terms which are of leading order near the boundary have 
been kept. Others have been neglected which, though smaller near the boundary, are of the 

42 



same size as some of the terms kept when integrated over the full area of the billiard. In 
particular, as it is, Eq. fl39|) gives the normalization of the wavefunction only up to C/kpA 
corrections. If however pw{E){l ± C/kpA) is used rather than pw{E) for the smooth part 
of the density of states, Eq. ([35D is replaced by 



1 ^ ^ ^ 1 ± jQ(2kpx) T 

" -Im Gosclr, r, E)^j^ —posc{E) 



(CI) 

Integration over the r.h.s gives precisely zero with each pair of terms from the beginning in 
order respectively canceling each other. Thus, the eigenstate normalization is unity from 
the leading constant term on the l.h.s. 
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